Energies of ions in water and nanopores within Density Functional Theory 
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Accurate calculations of electrostatic potentials and treatment of substrate polarizability are crit- 
ical for predicting the permeation of ions inside water-filled nanopores. The ab initio molecular 
dynamics method (AIMD), based on Density Functional Theory (DFT), accounts for the polariz- 
ability of materials, water, and solutes, and it should be the method of choice for predicting accurate 
electrostatic energies of ions. In practice, DFT coupled with the use of periodic boundary conditions 
in a charged system leads to large energy shifts. Results obtained using different DFT packages may 
vary because of the way pseudopotentials and long-range electrostatics are implemented. Using max- 
imally localized Wannier functions, we apply robust corrections that yield relatively unambiguous 
ion energies in select molecular and aqueous systems and inside carbon nanotubes. Large binding 
energies are predicted for ions in metallic carbon nanotube arrays, while with consistent definitions 
Na + and CI - energies are found to exhibit asymmetries comparable with those computed using 
non-polarizable water force fields. 
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I. INTRODUCTION 

As computational power increases, Density Functional 
Theory (DFT)-based aqueous phase ab initio molecu- 
lar dynamics (AIMD) simulations are becoming more 
widely used. DFT accounts for the polarizability of wa- 
ter molecules and can predict ion hydration structures 
significantly different from those computed using non- 
polarizable force fields. It also accounts for the polariz- 
ability of synthetic nanopores like carbon nanotubesii 2 " 
and the surface charge/dipole distribution- induced elec- 
trostatic potential inside inorganic nanopores, 3, which 
may strongly affect ion permeation, rejection, partition- 
ing, and exchange in these water-filled systems. The abil- 
ity to use AIMD to compare the intrinsic free energies 
of ions in bulk water and inside water-filled nanopores 
may revolutionize our understanding and modeling^ of 
ion transport in nanopores. The predictions can poten- 
tially be tested using recently synthesized membranes 
made of vertically aligned carbon nanotube^ or func- 
tionalized silica nanopores. 9 But up to now, when charge 
states of simulation cells are varied in the liquid phase, 
as when adding an ion to a simulation cell, the anoma- 
lous average, absolute, "intrinsic" electrostatic potential 
^10,11,12,13 nag n i nc i erec } direct AIMD determination of 
energies j 14 ' 15 i 16 

In solid state DFT calculations, when Fourier trans- 
form or Ewald techniques are used to compute long range 
electrostatic interactions in periodically replicated unit 
cells, it has been widely accepted that <f> is in general 
ambiguous i 1 ^ 11 ' 12 ! 13 So is the energy incurred q<p when 
adding an ion of charge q to a charge neutral simulation 
cell. Thus, successful estimates of the energy difference 
between charged species in vacuum and those embedded 
in a solid, such as the work function of electrons^ or 
in electrochemistrical applications^ have generally in- 
volved computing the difference in electrostatic poten- 



tials across an explicit interface in a single simulation cell. 
Since this difference depends on the interface, which is 
not always easily specified in complex systems, methods 
to reference the intrinsic electrostatic potential in solids 
to some known values, like the vacuum or Fermi level, 
have generated considerable interest ed 9 ' 20 ! 21 

In contrast, in physical chemistry, techniques for using 
classical force fields to compute the intrinsic free energy 
of ions in liquid water have been well established and 
widely applied^^^^^^^^i^ 2 . Here the elec _ 

trostatic component of the ion hydration free energy is 
computed using thermodynamic integration (TI) or free 
energy perturbation (FEP) techniques^ 3 - It is added to 
the short-ranged repulsion and van der Waal's contri- 
butions, which are often represented by Lennard- Jones 
terms. 

Formally, the predicted ion hydration free energies can 
be compared with thermodynamic estimates of the abso- 
lute hydration free energies of ions which are based on sol- 
vation and/or gas phase ion cluster dataj 34 i 35 To this end, 
the pure solvent contribution to the water-vapor electro- 
static "surface potential" 36 i 37 i 38 i 39 Acj) also needs to be 
computed. Separating the calculation into intrinsic and 
interfacial contributions is advantageous because electro- 
static effects are additive, and the pure solvent A(j> de- 
pends only on the water force field usedi 38 i 39 Many widely 
used water models yield Atfi on the order of ~ 1 kcal/mol 
per proton charge. 38 ' 40 ' 41 At water-material interfaces, 
A<p will include the substrate contribution. At vacuum- 
material interfaces, the substrate A<j) directly contributes 
to the work function and can be tested against photoelec- 
tric measurements! 17 i 42 As will be shown, incorporating 
parts of the substrate or pure solvent A<p is crucial for a 
consistent definition of intrinsic DFT ion energies. This 
is our main interest in the Aip term. 

Note that there is also an ionic contribution to aque- 
ous phase A(j) which does not vanish even at infinite 




FIG. 1: 

Samples of the model systems examined in this work, (a) 

CH 4 molecular crystal; (b) (6,6) SWNT fragment; (c) 
infinitely long (18,18) SWNT; (d) liquid water. The Na+ 
(CP) ion is depicted in blue (green). 

ion dilution^ 3 - Further, it has been suggested A</> at the 
water- vapor interface cannot be measured ! 36 ! 44 Neverthe- 
less, the combination of intrinsic ion free energy differ- 
ences in bulk and confined water and the interfacial A<f> 
may force ions into or reject them from narrow nanopore 
entrances and lead to novel ion blockage behavior at fi- 
nite electrolyte concentration. These effects will likely 
affect ionic currents through nanoporous membranes. 

In this work, we will use DFT to compute the intrinsic 
energy of ions in some condensed phase systems with- 
out explicit reference to an interface, just as has been 
done for ions in water using classical force fields. This re- 
quires regularizing the DFT electrostatic potential, which 
will be performed in two different environments: dipole- 
free molecular solids and metallic carbon nanotubes, and 
bulk liquid water. See Fig. [TJ The philosophies differ, 
but technically they both involve going beyond the well- 
known charge-image and charge-neutralizing background 
correctio n 19 i 20 i 23 and considering the "spherical second 
moment" electrostatic potential contributions of the con- 
densed phase materiah 13 ' 19 ! 37 Unlike pure Ewald repre- 
sentation of electrostatics, our corrections involve mak- 
ing physical choices, such as preserving the integrity of 
molecules and nanotubes. To facilitate this correction in 
liquid water, the electronic charge density is decomposed 
on to molecular centers using maximally localized Wan- 
nier functions ! 45 ' 46 To compare with DFT ion hydration 
results, we will focus exclusively on classical force field 
simulations that also apply Ewald sums. However, issues 



arising from truncating the coulomb interactio n 2 ' 7 i 8 
may also be relevant to DFT calculations that do not 
use Ewald sums. 

This paper is organized as follows. Section 2 briefly 
describes the computational methods. It reviews how 
(j) may be referenced to the vacuum or a known inter- 
face value, and emphasizes and reconciles the approaches 
taken in the solid state physics and physical chemistry 
communities. The results are presented in Sec. 3, and 
Sec. 4 concludes the paper with further discussions. An 
appendix describes the construction of maximally local- 
ized Wannier functions within the projected augmented 
wave (PAW) framework. 49 ! 50 

II. METHOD AND THEORY 
A. Systems without dipole moments 

We consider 

-^ion — -^ion+substrate -^barc ion -^substrate- (1) 

In periodic boundary condition DFT calculations, £"i on 
contains the generally unknown term qcj) , where q is 
the charge of the ion. This ambiguity arises from the 
inherently long-ranged nature of the Coulomb interac- 
tion V(r) = l/|r|, whose Fourier representation V(k) = 
47r/|k| 2 is undefined at k = 0. Ewald-sum calculations 
almost universally and arbitarily set V(k = 0) ex <f> to 
zero. 

In this work, we assume that any ion-containing pe- 
riodically replicated simulation cell has a neutralizing 
background, in effect removing the monopole. At zero 
temperature, the correction to E- lm due to the back- 
ground and the image charge interactions has been well 
documented in the DFT literaturei^ 9 - It is 

£cwaid = q 2 a/{2Le) (2) 

in atomic units, where a is the Madelung constant, L 
is the length of the cubic unit cell, and e = or e a 
depending on the application. When computing ion hy- 
dration free energy in liquids using classical force fields, 
a similar correction yield system size independent hydra- 
tion free energies ^2^SLSS^3m. Equation [2] does 
not remove the electric field arising from the periodic 
images, which is an area of active research . 20 ! 51 How- 
ever, for the purpose of this work, we will assume the 
monopole issue is a solved problem, and will apply Eq. [2] 
and its generalization to non-cubic cells with the appro- 
priate e. Instead, we focus on the correction to the in- 
trinsic electrostatic potential which persists, and can be 
on the order of electron volts, even in the limit of a very 
large unit cell where -E e waid becomes vanishingly small. 
The origin of this shift in <p is discussed below and is 
further illustrated in Sec. IIII Al 
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In general, a simulation cell exhibits finite dipole (d), 
quadrupole (the second rank tensor Q), and "spherical 
second" (SSM, Q) moments: 

d = y"rfr(r-R c )p(r); (3) 

Q = (1/2) J dr(r - R c ) : (r - R c )p(r); (4) 

Q = (l/3)TrQ. (5) 

Here R c is an arbitarily chosen center of the unit cell, 
p(r) is the total charge density, and the integration over 
r in the unit cell includes both electronic and nuclear 
degrees of freedom, 



J dvp{v) -+ J 



dr c p e (r e ) 
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,p8(r - R m ,a), 



where (3 denotes the nuclear site on molecule m, and e 
stands for electrons. 

Unit cells with non-zero d and or Q elements exhibit 
(fio that depend on the surface terminations. Consider a 
macroscopic but finite object built from a periodically 
replicated unit cell, terminated in interfaces with the 
vacuumJ^ If the lowest multiple moment of the the unit 
cell is a monopole, dipole, or quadrupole, the center of 
the object ("crystal") experiences an electrostatic poten- 
tial of the form J drr 2 f(r) exerted by the rest of the 
object. Here \f(r)\ ~ 1/r, 1/r 2 , or 1/r 3 respectively, 
and the vacuum region outside the object, r > R, is as- 
signed an electrostatic potential of zero. The integral is 
not uniformly convergent at large R, and thus <j) in unit 
cells at the center of the object depends on the shape of 
the object in addition to the nature of the interface. 13 
In particular, if |d| ^ 0, Q will depend on the choice of 
R c , and there will be an d-dependent electric field across 
each unit celli&lLlZdl 

Ewald methods, when applied with the "tin- foil" 
boundary conditions (almost universally the case in DFT 
calculations), effectively impose an infinite eoo dielectric 
continuum boundary that eliminates d- and Q-induced 
electric fields and potentials. As mentioned above, it 
sets <p — V(k = 0) = independent of R c . 13 To report 
E- lon in the absence of an interface, the "intrinsic" non- 
zero q<j) , which depends on physical choices of system 
boundaries, needs to be added to the Ewald-computcd 
ion energy. Such a regularization of DFT electrostatic 
potential within Ewald sum calculations has been carried 
out in Ref. [l^ for an isolated ion in the gas phase. The 
importance of this correction has not been sufficiently 
emphasized in solid or liquid systems. 

We first consider special cases where d vanishes or 
can be chosen to vanish, and obvious physical choices 
of simulation cell boundaries can be defined. Such sys- 
tems include ions inside a crystal made up of molecules 



without permanent dipole moments, and ions embed- 
ded in carbon nanotubes. Consider the spherical sec- 
ond moment (SSM) contribution in these systems. This 
term docs not appear in real space multipole expan- 
sions of the electrostatic potential exerted by an arbitary 
charge distribution^ 3 - But when Ewald sums arbitarily 
set (j) a = 0, an SSM-dependent constant shift arises. This 
shift is undone by adding ; 13 ' 19 



issM(Rc) = -(2tt/30) J dvp{v){v - R c ) 2 



(6) 



where is the unit cell volume. If p(r) can be de- 
composed into molecular contributions labeled by "m" 
and centered on molecular centers R m , such that p(r) — 
J2 m Pm{r), then 

0ssm(R c ) = -(2tt/30)^ f drp m (r)(r ~ R m ) 2 



-2 J dr(R c - R m ) • (r - K m )p m (r) 



(7) 



where we have used the fact that each molecule has no net 
charge. The position of the ion, if present, is assumed to 
coincide with R c , which permits the use of Eq. [7] without 
adding monopole terms. If each molecule also possesses a 
zero dipole moment, the second term vanishes, and Eq.[7| 
simplifies to 
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(8) 



In this case, </>ssm no longer depends on R c . A system 
of identical, non-interacting molecules without dipole 
(d m = 0) or quadrupole moments, ([Q m — l/3TrQ m ]jj = 
0,* + j), Qm = (1/2) J drp m (r)(r - R m ) : (r - R m ), 
would then yield a </>ssm proportional to the particle den- 
sity and the molecular-centered term J drp m (r)\r— R m | 2 . 
The absence of d m means that Q m is also independent 
of R m . 

Several points are worthy of note. (1) Unlike pure 
Ewald summation, applying SSM corrections amounts to 
making physical choices of system boundaries. The elec- 
tron density p e (r) is not uniquely assigned to a simulation 
cell, but can be varied via "spreading transformations."— 
For example, one can split the p e (r) of a molecule which 
straddles the unit cell boundary, creating two or more 
charge fragments; the 0ssm correction, now no longer 
decomposable into intact molecular contributions, will 
depend on how this splitting is performed. Physically 
meaningful interfaces consist of intact molecules without 
splitting of molecular orbitals across unit cell boundaries. 
When modeling ions in nanotubes, the choice of bound- 
aries will clearly affect the qfissM corrections to ion bind- 
ing energies inside the pores. Here the simulation cells 
should preserve the p e (r) of intact SWNT in the lateral 
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directions. Splitting of the continuous p e {r) along the 
axial direction is necessary, but intuitive choices can still 
be made to demarcate p e {f) within each unit cell of an 
infinitely long SWNT, guided by symmetry considera- 
tions. In covalently bonded solids like silicon, the choice 
of the unit cell is less obvious. (2) Note that carbon 
nanotubes are generally chemically modified at their two 
open ends, leading to surface contributions that may need 
to be added to the computed intrinisic electrostatic po- 
tential when predicting ion transport through synthetic 
nanoporous membranes Ionic contributions will also 
be present^ 3 . (3) </>ssm depends on whether core elec- 
trons or pseudopotentials are used via Eq. [S] Correct- 
ing <j) with (/<ssm is thus critical for comparing intrin- 
sic ion energies computed using different DFT codes and 
methods. (4) For unit cells without dipole moments, but 
with finite quadrupole moments, [Q — l/3TrQ]ii ^ 0, 4> Q 
still mathematically depends on the shape of the macro- 
scopic system, but not on unphysical considerations such 
as whether core electrons are used. (5) For isolated ionic 
species in vacuum, q<f>sSM vanishes when f2 — > oo. In 
condensed phase systems, where finite electronic and nu- 
clear charge densities persist in all space, q<t>ssM becomes 
f2-independent for large simulation cells, and is roughly 
proportional to the atomic density. As will be seen, it can 
result in a several volt correction to the DFT-computed 
<j) - (6) As will be discussed, classical force fields also 
exhibit <?</>ssm corrections to the electrostatic potential, 
albeit smaller ones that contain no ambiguity associated 
with "core electrons" contributions. 



B. Ions in liquid water 

We also consider the energy of ions in bulk liquid water, 
where Eq. |T| is averaged over the AIMD trajectory. Here 
E- lon is not the ion hydration enthalpy. It does not ac- 
count for the solvent reorganization penalty arising from 
molecular motion, and is instead related to the value of 
the integrand at the end point of the thermodynamic in- 
tegration formula where the ion charge q\ is gradually 
increased from zero to q£2- See Sec. IIII CI See Sec. Ill 
C. However, regularizing cf> is still crucial for correcting 
Eion(qx) along the entire QA-integration path. As such, 
our work lays the foundation for future AIMD thermo- 
dynamic integration calculations. 

In general, a simulation cell containing an instanta- 
neous configuration of liquid water has non-vanishing 
dipole and quadrupole moments. Each water molecule 
also has finite d m and Q m elements. The finite d m means 
that Q m now depends on the choice of the molecular 
center^ Despite this, following formulations developed 
for classical water force fields (see below) , <?</>ssm (Eq. [8]) 
is still applicable because neither the vectors (R c — R m ) 
nor the molecular axes (implicit in J dr) has preferred 
orientations, and the second term in Eq. [7] averages to 



zero in liquid environments. We will collapse the elec- 
tron density into molecular contributions via the max- 
imally localized Wannier function approach ) 45 i 46 taking 
care to preserve the integrity of H2O molecules and not 
cleaving orbitals across unit cell boundaries. 

Formally, comparing predicted absolute ion free ener- 
gies with thermodynamic estimates based on solvation 
and/or gas phase ion cluster dat a 34 ' 35 requires adding 
the pure solvent contribution to A<j> at the water-vapor 
interface (henceforth simply referred to as A<f> in this sec- 
tion). The A0 formulation is well established in the clas- 
sical force field literature ! 36 ' 37 

Acj, = /47r[-(l/3)ApTrQ ro + ^ ' p d (z)]\ (9) 
= A(j)Q + A<j) d 
Pd{z) = (^J drp m (r)5[z-(R m ) z }(d m ) z y (10) 

Here z is the direction perpendicular to the interface, Ap 
is the difference between water and vapor densities, and 
Pd(z) is the dipole density profile. d m and Q m are the 
dipole vector and quadrupole moment tensor for water 
molecule m. They are in fact independent of m because 
of liquid phase self-averaging ("(...)"). Equation would 
involve a sum over molecular species if more than one 
were present. Its derivation makes use of the lack of 
preferred molecular orientations, mentioned above. Al- 
though Refs. |H and H3, which preceed Ref. [H, do not 
use the term "second spherical moment," the first term 
on the right side of Eq. [5] is clearly equivalent to 4>ssm 
defined in Eq. [5] 

Thus, the intrinsic ion hydration free energy computed 
using classical force fields is typically defined without the 
q<fisSM term, which is instead included as part of the pure 
solvent surface potential A(/>j 55 i 56 i 57 The numerical values 
of pd{z) and Q m depend on the choice of the molecular 
center R m , but E- mn +qA(j) (or more precisely the free 
energy AGi on + qA<p) is independent of such choices.— 

If we compute qA(j> using an AIMD simulation of the 
water-vapor interface^, and add it to an appropriately 
defined AIMD intrinsic ion hydration free energy, the 
result will be unambiguous. However, using AIMD to 
compute A<f> is not yet feasible because (a) large numeri- 
cal noises are associated with the liquid- vapor interfacial 
fluctuations, and overly long and costly AIMD trajec- 
tories would be required to yield good statistics for A(j>. 
(b) Widely used DFT functionals do not necessarily yield 
1.0 g/cc water density 58 ' 59 rendering the use of Eq. 
which depends on the difference between bulk water and 
vapor densities Ap, problematic. 

At present then, we believe AIMD predictions of ion 
free energies are most useful for comparison with and 
validation of classical force fields. This comparison re- 
quires that the q(f>ssM correction be consistently included 
as part of the intrinsic ion energy, and not as part of A</>, 
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for both AIMD and classical force field simulations. In 
fact, if they are both perfectly accurate and consistently 
defined (using the same molecular centers R m ), DFT and 
classical force field calculations of A(j>d (Eq. [T0|) , and not 
just the intrinsic ion energies, should agree with each 
other. This is because pd{z) only depends on the spatial 
and orientational distributions of water molecules and 
their instantaneous dipole moments at the water-vapor 
interface. But Acf> will not be the focus of this work; our 
interest in this interfacial quantity is mainly to establish 
a consistent definition of E lon such that DFT, force field 
predictions, and thermodynamic dat a 34 ' 35 can be com- 
pared. 

Recall that E lon neglects solvent reorganization cost 
arising from molecular motion - the dominant contri- 
bution to hydration in liquid water. Using Ewald-sums 
and non-polarizable classical force fields, this definition 
of .Eion (or "q<fi" in the notation of Ref. 23) becomes inde- 
pendent of the simulation cell size when it is corrected by 
twice the value of Eq.[2j 23 We will apply this correction to 
AIMD-predicted E lon to compare with force field results. 
Here the DFT electron density already "reorganizes" to 
accommodate the ion and introduces a slight ambiguity. 
The Ref. [23| correction will be shown to be reasonable by 
testing -Eion convergence with system size. In contrast, 
the absolute ion hydration free energy (AG), a physical 
quantity, contains both electronic and nuclear solvent re- 
organization. It should be unambiguously cell-size inde- 
pendent when Eq. [5] is subtracted once^ As mentioned 
above, for ions inside carbon nanotubes, electronic relax- 
ation dominates; we will freeze nuclear degrees of freedom 
and show that Eq. [5] suffices for convergence. 

To recapitulate, we consider two types of systems: ions 
in dipole-free unit cells and in liquid water, both perti- 
nent to ion penetration into water-filled nanopores. Dif- 
ferent philosophies, traditionally associated with solid 
and liquid phases, are involved, but the </>ssm correc- 
tions (Eqs. [5] and [5J are crucial in both cases. For a 
system of dipole- and quadrupole-free molecules, the two 
viewpoints coincide. 4>Q (Eq. [9]) becomes independent 
of the choice of molecular centers, Pd{z) vanishes iden- 
tically, and Eq. [5] reduces to Eq. [5J In this case, Eq. [5] 
or Eq. [8j directly references Ei on to ^Vacuum = through 
a material- vacuum interface that has no dipole contribu- 
tion to A(j) (Fig. UK). 



C. VASP calculations 

We apply the Vienna Atomistic Simulation Package 
(VASP)^o and the standard VASP suite of ultrasoft 
(US) 61 and projector-augmented wave (PAW) pseudopo- 
tentials (PP)4&22, US PP are used in AIMD simulations 
while static calculations utilize PAW PP. The Perdew- 
Wang (PW91) and Perdew-Burke-Ernzerhof (PBE) ex- 
change correlation functional a 62 i 63 are applied for ions in 



bulk liquid water and ions inside molecular solid/carbon 
nanotubes, respectively. These functionals generally 
yield similar predictions. The Na + PP does not in- 
clude pseudovalent p electrons. The AIMD time step is 
0.25 fs, the cubic simulation cells have sizes of 9.9874 A 
(12.509 A) for 32 (64) H 2 molecules plus one ion, all 
proton masses are replaced by deuterium masses, and 
the plane-wave energy cutoff is set at 400 eV. The con- 
vergence criterion at each of the Born Oppcnhcimcr 
time steps is 10 -5 eV. The total energy is conserved 
to within 1 K/ps (2 K/ps) for the CU (Na+) simu- 
lation cells, respectively. An elevated temperature of 
T=375 K is enforced using a Nose thermostat to alle- 
viate water-overstructuring problems when applying the 
PW91 functional for AIMD simulations.^ 4 .! 6 ^ The tra- 
jectory lengths used, not counting pre-equilibration runs, 
are at least 20 ps. In one case (CI - in 64 H2O box), we 
actually use two different initial configurations obtained 
from classical force field molecular dynamics simulations, 
conduct two 10 ps AIMD trajectories, and combine the 
results. The mean E- lon of the two AIMD segments are 
within 0.034 eV of each other. In the aqueous phase cal- 
culations, i?i on (Eq. [J) is sampled every 0.1 ps. T-point 
Brillouin zone sampling is applied in most cases, except 
that a 1x1x10 k Monkhorst-Pack grid is also used in 
convergence tests for infinitely long metallic carbon nan- 
otubes. 

The implementation of the maximally localized Wan- 
nier functions in VASP is limited to a T-point only 
sampling of the Brillouin zone. Details concerning the 
construction of maximally localized Wannier functions 
within VASP are provided in the Appendix. 

To apply Eqs. II] and [H the VASP subroutine pot.F 
is slightly modified so that the electrostatic potential at 
G = 0, k = reflects the total number of electrons (not 
the negative of the total pseudo- nuclear charge). 



III. RESULTS 
A. Ions in molecular crystals 

We first consider two test cases of hypothetical simple 
cubic crystals made up of Ar atoms and CH4 molecules. 
A lattice constant of 5 A ensures that the molecules are 
well-separated. We construct a 3x3x3 supercell, create 
a vacancy at the supercell center, and insert a Na + or 
Cl~ ion there. The Ewald corrections (Eq. [5]) to elim- 
inate the effects of the neutralizing background cancel 
if we assume e = 1 for both the isolated ion and the 
ion inside the molecular crystal vacancy^ Centering the 
ion in this supercell ( "unit cell" ) leads to zero net dipole 
and quadrupole moments. This physical choice of cell 
boundary allows unambiguous computation of q4>sSM us- 
ing Eq.[5] (i.e., without decomposing p(r) into molecular 
contributions). Furthermore, for a macroscopic crystal 
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system 


cell size 


Na+ 


cr 


CH 4 (s) 


15.0x15.0x15.0 


-0.22 (-1.33) 


-0.33 (+1.00) 


Ar(s) 


15.0x15.0x15.0 


-0.12 (-1.66) 


-0.20 (+1.48) 


C48H24 


15.0x15.0x15.0 


-1.11 (-2.48) 


-0.10 (+1.51) 


C48H24 


20.0x20.0x20.0 


-1.18 (-1.76) 


-0.15 (+0.53) 


(6,6) 


14.8x14.8x14.8 


-2.85 (-6.29) 


-1.77 (+2.14) 


(6,6) 


14.8x14.8x22.2 


-2.86 (-6.26) 


-1.80 (+2.16) 


(18,18) 


17.0x17.0x17.0 


-2.20* (-7.04) 


-1.09 t (+4.11) 



TABLE I: Ei on inside molecular crystals and SWNT's, in units 
of eV. Simulation cell sizes are in A 3 . The values in brackets 
are computed without the </0ssm correction. *K + rather than 
Na + is used in the (18,18) SWNT; ' overestimated due to a 
small amount of electron transfer (see text). 



made up of this unit cell, any crystal-vacuum interface 
will not contribute to <fi , A(/>=0, and our predicted Ei on 
are unambiguously referenced to the vacuum value = 0. 
No geometric relaxation of the lattice atoms is allowed 

for the isolated 



are ion 



after inserting Na + or Cl _ . E\ 

CI" in a (15 A) 3 unit cell also exhibits a small qifissM = 
—0.1 eV using pseudopotential containing only 3s and 3p 
valence electrons. This is subtracted from the CI "-plus- 
substrate g^sSM term. With these corrections, we obtain 
E ion =-0.22 (-0.33) eV for Na+ (CI - ) in the CH 4 solid, 
and E- lon =-0.12 (-0.20) eV in Ar solid. These values are 
reasonable. With these large band gap atoms/molecules, 
we expect small interactions between the ions and the 
molecular dipoles they induce in the substrate. 

More significantly, adding g</>sSM leads to comparable 
.Eion for both ions. This is expected because the bind- 
ing energy should be roughly proportional to q 2 . In the 
absence of q4>sSM corrections, Na + and Cl _ ions exhibit 
Ei on ~ ±1 eV in our model CH 4 and Ar solids (see Ta- 
ble [Q). The large asymmetry between cations and anions, 
and the predicted repulsion of Cl~ , mean that these un- 
corrected DFT energies are clearly unphysical. Table U 
implies that ^ssm > in the DFT calculations. This is 
because nuclei are extremely localized in space whereas 
electrons are more diffuse. The true electrostatic poten- 
tial <j) averaged over the unit cell includes contributions 
inside atomic cores where <f> > 0. Thus <p is positive def- 
inite inside the crystal, not zero as Ewald sums typically 
mandate. The position independent </>ssm term correctly 
accounts for this. Similar corrections have been applied 
to non-polarizable Ne solid and Ar fluid models using 
analytical atomic charge distributions ! 21 ' 56 

The 0ssm correction is a consideration whenever clas- 
sical force fields with distributed charges and Ewald 
sums are used to compute intrinsic hydration free ener- 
gies of ions i 68 ' 69 Even without distributed charges, clas- 
sical force field descriptions of molecules still exhibit 
such corrections. If the C and H atom sites in our 
CH4 solid carry fixed point charges of — 0.4e and +0.1e, 
<pSSM = —0.11 volt. Adding q(f>ssM will then yield the 



expected physical result, namely, that the ion energies 
in the vacancy site are almost zero when computed us- 
ing non-polarizable force fields. Recall that </>ssm can 
alternatively be defined as part of the interface surface 
potential (Eq. [5]) in force field calculations. 



B. Ions in carbon nanotube models 

Unlike the above molecular crystal test cases, the 
interactions of ions inside single wall carbon nan- 
otubes (SWNT) are of significant technological inter- 
est. Narrow pore (~ 8 A diameter) SWNT arrays 
have been proposed as osmotic membranes^ and axially 
aligned SWNT's with larger pore diameters have recently 
been synthesized However, most molecular dynamics 
simulations of water-filled SWNT^M ignore SWNT 
polarizability^ which is infinite for metallic tubes, and 
large polarizability-induced effects are expected. 

Figure |T|b depicts a proton-terminated (6,6) SWNT 
fragment C48H24 situated in the center of a 15 x 15 x 15 A 3 
cubic unit cell. This finite-sized "molecule" has a finite 
band gap and no dipole moments. We place a Na + or 
a Cl~ ion at the center of the C48H24 interior. Ewald 
corrections (Eq. [5]) used to eliminate monopole contribu- 
tions are again assumed to cancel for the bare ions and 
the ion-SWNT fragment complexes (see below). The net 
Eion, computed without relaxing SWNT atoms, are pre- 
dicted to be -1.11 eV and -0.10 eV for Na+ and Cl~, 
respectively (Table [}. Without qfissM (Eq. O, they are 
erroneously predicted to be -2.6 eV and +1.6 eV. In- 
dependent verification of this ion-position independent 
correction is obtained by moving the Na + from 
the tube interior to an axial position halfway between 
the tube openings. The energy difference between these 
configurations is 1.17 eV. Thus, E ion ~ —1.11 eV indeed 
reflects the binding energy of Na + in C48H24. We em- 
phasize that, using non-polarizable force fields* 3 -! 5 -' 6 - the 
binding energy would be purely dispersive and would at 
most be a small fraction of an eV. 

The binding energy should be independent of the sim- 
ulation cell size for this isolated SWNT fragment. With 
a larger, L = 20 A, cell the q4>ssM corrected Na + (Cl~) 
binding energy becomes 1.18 eV (0.15 eV). The small 
increase compared to L = 15 A may be due to the po- 
larizability of C48H24 in the finite simulation cell, which 
may make eoo slightly larger than unity. As L — > 00, 
both q4>ssM and (e^ — 1) vanish as 1/L 3 . Nevertheless, 
|<70SSm| remains substantial, ~ 0.6 eV, for the L = 20 A 
cell (Table [XJ) . It is much larger than that derived for 
isolated monoatomic ionsi^ because of the sheer number 
of electrons present in C48H24. 

Next, we consider Na + /Cl~ in infinitely long (6,6) and 
(18,18) SWNT arrays. These models mimic micron-thick 
carbon nanotube sieves which have been synthesized. 8 
The SWNT's chosen in the present study are metallic and 
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exhibit — > oo in the axial direction. These nanotubes, 
with cubic unit cell sizes chosen to reflect the SWNT 
periodicity (14.8x14.8x14.8 A 3 and 17.0x 17.0x 17.0 A 3 , 
respectively), are not finite molecules. However, by cen- 
tering the ion inside the SWNT and the SWNT within 
the simulation cell, and truncating the electron density 
at the z-direction of the grid boundary — taking care to 
evenly distribute p{r) at the z-edges — we have made a 
physical choice of p(r) in the unit cell that exhibits no 
dipole moment and permits qcf>sSM corrections. The in- 
trinsic energies of ions as they are displaced from the pore 
center can then be referenced to the pore center value. 

Equation [2] corrects the energies of the isolated ions. 
Ions inside the SWNT are assumed to have E ewa \ c i = 
because of the large eoo of the SWNT array. Table Q] lists 
the resulting E- mn for Na + and CI - . Na + and CI - ex- 
hibit large binding energies of 2.9 and 1.8 eV respectively 
in the narrow (6,6) metallic SWNT arrays. Increasing 
the Brillouin zone sampling from T-point to a 1x1x10 
Monkhorst-Pack grid leads to changes in i?i on that are 
less than 0.02 eV. 

Unlike the C48H24 fragment, the model SWNT arrays 
are meant to be periodically replicated in the lateral di- 
rections. We consider system size dependences only by 
increasing the axial dimension of the (6,6) simulation cell 
to 22.2 A. We find that E lon are conserved to within 
0.03 eV after adding g^sSM, which are on the order of 
several eV (Table This suggests that our somewhat 
arbitary truncation of electronic density at the z-edges 
nevertheless leads to robust (?0ssm corrections. Unlike 
the Ar and CH4 crystals, where A(j> = 0, E- lon cannot 
yet be referenced to vacuum values because we have not 
specified the physical interface that will determine A<f>. 
However, the two qA(f> should cancel when we consider 
the sum of Ei on for the oppositely charge Na + and Cl~ . 
Thus, regardless of the surface terminations, we pre- 
dict substantial binding energies for ions inside metallic 
SWNT. 

No charge transfer occurs between Na + or Cl~ and the 
(6,6) SWNT array. DFT/PBE does predict that Na+ ac- 
cepts a small fraction of an electron from the (18,18) ar- 
ray. As a result, we switch to the less electronegative K + , 
whose lowest unoccupied orbital is found to reside above 
the (18,18) Fermi level. Thus, K + is representative of a 
monovalent cation well separated from the SWNT pore 
surface. It exhibits a still considerable 2.20 eV binding 
energy inside the wider (18,18) array. CI - experiences a 
1.09 eV binding energy. This slightly overestimates the 
interaction of this SWNT with a monvalcnt point an- 
ion because DFT/PBE predicts a small charge transfer 
(~ 0.3 electron) to the SWNT. 

Studying ions inside SWNT is pertinent to permeation 
of electrolytes, where the presence of water strongly sta- 
blizes Na + and Cl + and helps prevent electron transfer to 
or from the nanotubes. How the strongly confined water 
will screen the induced electrostatic interaction between 




Hps) 

FIG. 2: 

-Eion as functions of time for (a) Na + and (b) CI - in a 
simulation cell with 32 H2O. There is no significant drift 
over these short 20 ps trajectories. 

ions and SWNT will be examined in the future, but the 
SWNT polarizability clearly will have a large impact on 
ion permeation, especially in the narrow (6,6) SWNT 
that accommodates only a single file water inside^^ 



C. Ions in water 

In this section, we compute the mean E{ on for Na + 
and CI - in liquid water. We apply the maximally local- 
ized Wannier function metho d 45 ' 46 to decompose electron 
density on each H2O. This method correctly locates four 
Wannier orbital centered within ~ 1.5 A of each water 
oxygen site, corresponding to the two O-H covalcnt bonds 
and two lone pairs. We choose O as the molecular center 
in accordance with classical force field conventions^ and 
add Ewald corrections (twice that of Eq. [21 as discussed 
above) and q4>sSM to the instantaneous E lon depicted in 
Fig. O -Eion exhibits little drift over these 20 ps trajec- 
tories. 

Before we examine E{ on further, the robustness of using 
Wannier functions to compute q<f>ssM corrections (Eq. [5]) 
needs to be documented. Over a sample of 14 configura- 
tions, evenly distributed over the Na+ AIMD trajectory, 
we find that <?</>ssm experiences only ~ 0.01 eV variation 
at each time step, amounting to a negligable standard 
deviation of 0.0028 eV. 

Our maximally localized Wannier functions, by con- 
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0.118 




0.01 



FIG. 3: 

<^>ssm (Eq. El in units of volt per H2O per (10 A) 3 ) as the 
simulation cell size L varies. Crosses: 0ssm for an isolated 
H2O, using maximally localized Wannier function. Solid 
line: linear fit to scaling behavior. Circle: <^ssm (Eq. EJ 
computed using p(r) directly (without using Wannier 
functions). Diamond and dashed line: </>ssm for 32 H2O 
(L = 9.9874 A) and 256 H 2 (L = 19.9748 A)3 



struction, depend somewhat on the simulation cell size4£ 
So do quantities we use to compute <^ssm within the Wan- 
nier function framework. For a single water molecule in 
a cubic cell of length L, the entire charge density be- 
longs to the one H2O, and the Wannier function-derived 
Eq. [6] can be compared with the value obtained from 
numerical integration of p{r). Figure [3] shows that the 
Wannier-derived 0ssm convergence oc L~ 2 ; </>ssm exhibits 
an empirical correction +0.623/L 2 volt- A 2 per H2O in a 
(10 A) 3 unit cell, or +623/L 5 volt-A 5 per H 2 0. The scal- 
ing is consistent with analysis of (r) and (r 2 ) described 
in the Appendix. (Recall </>ssm scales as l/il for iso- 
lated molecules.) Numerical integration over p(r) yields 
L-independent 4>ssm- Its value agrees to within 0.01 % 
of the Wannier 0ssm when the latter is extrapolated to 
L — > 00. This extrapolation increases g^sSM computed in 
a L ~ 10 A cubic cell by ~ 5.5 %. For water at 1.0 g/cc 
density (32 H2O molecules for this cell size), this increase 
amounts to 0.22 volt. 

Next, we show that the presence of multiple molecules 
does not substantially change the L-dependence. We ex- 
pand al = 9.9874 A box of 32 water molecules into a 
2x2x2 supercell, and compare the difference in <^>ssm for 
this water configuration. The L = 19.9758 A supercell 



ion 


iVw 


<£ion) 


cr 


32 


-6.57 (-2.60) 


cr 


64 


-6.27 (-2.29) 


crt 


256 


-7.39 (-8.23) 


Na+ 


32 


-7.85 (-11.54) 


Na+ 


64 


-7.70 (-11.36) 


Na +t 


256 


-10.01 (-9.17) 



TABLE II: Mean £ ion in eV, accumulated using the PW91 
functional with (without) the q<j>ssM correction. N w is the 
number of water molecules in the simulation cell. ^ SPC/E 
classical force field results (Table 4, Ref . I23I) . which are inde- 
pendent of simulation cell size. 



</>ssm exceeds the L = 9.9874 A value by 5 %. Accord- 
ing to the scaling curve in Fig. [31 the Wannier function 
approach yields a 4.2 % difference between L — 20 and 
L = 10 A for an isolated H2O. This suggests that extrap- 
olating q<pssM computed for a L ~ 10 A 3 box of liquid 
water to L — > 00 at 1.0 g/cc density, using the L-scaling 
for a single H2O, should only lead to a small systematic 
error of about (5-4.2)/5 x 0.22 volt, or about 0.035 volt. 

Having shown that g^ssM corrections generated using 
maximally Wannier functions are robust, and that the 
errors are well-controlled, we return to E- lcm in liquid wa- 
ter (Fig. [2]). E ion values averaged over at least 20 ps are 
listed in Table II. q^ssm is computed using Eq. 8 in the 
presence of the ion, with (?0ssm for the isolated Cl~ ac- 
counted for. Removing the ion before computing q<pssM 
changes this term by less than 0.01 eV. Overall, (j^ssm is 
large — on the order of several eV — for the water-filled 
simulation cell. Na + and Cl~ exhibit a 1.28 eV asym- 
metry in £ion- To properly compare this with force field 
results, we also apply the g^sSM correction to SPC/E 
water model used in Reb^ 3 - with oxygen as the molecu- 
lar center. Because of the SPC/E charge distribution, 
this increases {decreases) the magnitude of E{ on for Na + 
(Cl~). These force field based Na + and Cl~ E lon now 
exhibit a 2.39 eV difference, which is substantially larger 
than the DFT/PW91 result. The sum of the PW91 Cl" 
and Na + E lon: -14.42 eV, is also substantially smaller 
than the -17.40 eV predicted using the non-polarizablc 
SPC/E water model. 23 The comparisons depend on the 
force field used. 

-Eion (Eq. [1]) is the end point of thermodynamic inte- 
gration (TI) of ion hydration in systems without elec- 
tronic polarizability. With polarizable water force fields, 
it should be replaced by (dE- lon / ' d(Xqi on )) x=i , where gi on 
is the charge of a non-polarizable ion force field. AIMD 
treats water (and ion) polarizability, but this formula- 
tion is in general not applicable because <7i on is set by 
total number electrons in the system, and the highest 
occupied and lowest unoccupied orbitals may not gen- 
erally reside on the ion. Instead, a quantum mechan- 
ics/molecular mechanics (QM/MM) approach with the 
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ion treated classically, followed by free energy perturba- 
tion to replace the MM ion by a DFT ion, may be a 
preferred thermodynamic route to compute AIMD ion 
hydration free energies. 

Thus, while the ion hydration free energies reported 
in Ref. |23j are in reasonable agreement with experiments 
and the DFT results are smaller and might suggest that 
DFT/PW91 underestimates monovalent ion energetics in 
water, comparing the end-point integrand value of a ther- 
modynamic integration formula between non-polarizable 
force fields and fully polarizable, DFT-based AIMD 
methods may be misleading. Furthermore, thermody- 
namic theories that incorporate DFT ion-water clus- 
ter energies and dielectric continuum approaches have 
yielded good agreement with experiments^ 

Classical force field intrinsic Ei on exhibit O(l) kcal/mol 
cell size dependences after correcting for image and back- 
ground charge contributions, 2 -" However, there have been 
recent suggestions that AIMD simulations exhibit larger 
cell size dependences for charged species J£ Thus, we also 
consider a Na + /Cl~ plus 64 H2O AIMD simulation at 
1.0 g/cc H2O density for 20 ps. We find that E; on changes 
by +0.15 eV (+0.30 eV) for Na+ (CI - ). This amounts 
to a 3.5 (6.9) kcal/mol L-dependence. 

Are these L-dependences statistically significant? It 
is difficult to estimate the uncertainty in Ej on using our 
short, ~ 20 ps AIMD trajectories (Fig. [5]). Using the fact 
a 1 ns classical force field molecular dynamics simulation 
typically yields a standard deviation of ~ 0.5 kcal/mol 
in .Eion, we estimate that our 20 ps AIMD simulations 
have 3 kcal/mol statistical uncertainties. In this sense, 
the Eion difference between the 32- and 64-H2O cells is 
within the expected statistical errors. This test seems 
to suggest that our use of twice the value of Eq. H 2 x 
2.04 eV, to correct the monopole-image interaction leads 
to a reasonable system size convergence (see Sec. Ill BP . 
However, the AIMD solvent electronic polarizability also 
plays a role in modifying the system size dependence seen 
in Ref. [23]. Further studies on ion hydration free energies 
are needed to address this issue. 



IV. CONCLUSIONS 

In this paper, we apply DFT methods to compute the 
energies of ions in condensed phase systems. In general, 
due to the ambiguity in the average intrinsic electro- 
static potential <p arising from periodic boundary con- 
ditions and Fourier transforms/Ewald sums, this energy 
is ill-defined. It even depends on the number of elec- 
trons (valence vs. core) and computational protocol. For 
two types of systems which are pertinent to ion perme- 
ation in water-filled nanopores, we demonstrate that the 
intrinsic ion energy can nevertheless be reported. The 
key is to apply the "spherical second moment" correc- 
tion (^ssm) ) 13 ' 37 in addition to the well-known Ewald 



correction pertinent to a charged simulation cell with a 
uniform neutralizing background J^ 2 ^ A physical choice 
of unit cell boundary is involved in this correction. 

The first type of systems consists of molecules or nan- 
otubes whose unit cell exhibits no permanent dipole mo- 
ment. This enables relatively unambiguous estimates of 
and yields ion energies that can be referenced to 
the vacuum phase. Note that the ion energy will still 
mathematically depend on the shape of the crystal if the 
unit cell contains a non- vanishing quadrupole moment. 
Surface potential terms may also contribute if the ma- 
terial surface is chemically different from the bulk re- 
gion. Fortunately, interfacial effects are additive to in- 
trinsic ion energies. With this approach, and choosing 
cell boundaries that do not split up nanopore electronic 
densities in radial directions, we have computed the in- 
teraction between an ion and model molecular crystals 
and carbon nanotubes. The intrinsic binding energies 
for metallic (6,6) tubes are found to be 2.9 and 1.8 eV 
for Na + and Cl _ , respectively. Even the much wider 
(18,18) SWNT's exhibit binding energies of 2.2 and (a 
slightly overestimated) 1.1 eV for monovalent cations and 
anions, respectively. This suggests that the proposed 
use of vertically aligned SWNT sieves as ion rejection 
membranes^^ may need to be re-examined. 

The second type of systems consists of liquid-state sim- 
ulation cells with no obvious symmetry or other means 
to eliminate the net unit cell dipole moment. The phi- 
losophy is different here, and follow formulations derived 
for classical water force fields, but qifrsSM corrections are 
still useful because of orientational self-averaging in the 
liquid state. Here we apply maximally localized Wannier 
functions to decompose charge densities onto molecules. 
<^SSM will depend on the choice of molecular center, which 
needs to be consistently defined to compare with classical 
force field results. It also involves the physical choice that 
each water molecule is preserved intact with no splitting 
of electron density across cell boundaries. With consis- 
tent 0ssm protocols and using a 32-water cell, we find 
that Na + and Cl~ exhibit ion energies (i.e., hydration 
energies neglecting solvent reorganization) of -7.9 and - 
6.6 eV respectively. These are considerably smaller than 
classical force field ion energies of -9.0 and -7.1 eVM Us- 
ing a large simulation cell containing 64 H2O molecules 
leads to 0.15 and 0.30 eV reductions in the ion energies. 
These differences are within the statistical uncertainties 
of the simulations. Note that our "ion energy" in the 
aqueous phase neglects solvent reorganization and is not 
a measurable quantity. It also differs from the end-point 
of the thermodynamic integration formula for the hydra- 
tion free energy because of electronic polarizability issues. 
Nevertheless, with consistent ^ssm corrections, Na + and 
Cl _ ion energies become comparable, do not contain am- 
biguities associated core electrons in pseudopotentials, 
and their asymmetry is now comparable to force field 
predictions. The ability to correct for 4> a is important 
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for future ion hydration free energy calculations using 
AIMD in conjunction with thermodynamic integration 
methods. 

Looking ahead, it is tantalizing to contemplate how the 
methods described in this work can be applied to com- 
pute ion energies in crystalline systems and nanoporous 
solids.— 



to the basis of reciprocal lattice vectors. The weights ujj 
are defined as in Appendix A of Ref. [46|. 

In the PAW metho d 49 i 50 the one-electron wave func- 
tions \I/ n are derived from the pseudo-wave-functions ^ n 
by means of a linear transformation 

l*»> = l*n> + Ed^) - \4>kWk\^n), (A.5) 
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APPENDIX: MAXIMALLY LOCALIZED 
WANNIER FUNCTIONS. 

The maximally localized Wannier functions W n are de- 
fined as linear combinations of the iV occ selfconsistent 
Bloch wave functions \& n that make up the space of oc- 
cupied orbitals, 



\W n ) = £ EW|*n')> 
n' = l 



(A.l) 



where U is the N occ x N occ unitary transformation matrix 
that minimizes the spread functional, 



E[U] = {{W n \r 2 \W n ) - (W n \r\W n ) 2 ) , (A.2) 

n=l 

which provides a measure for the degree of spatial local- 
ization of the Wannier functions. (Since the implemen- 
tation of the maximally localized Wannier functions in 
VASP is limited to a T-point only sampling of the Bril- 
louin zone, we have dropped the usual k-point index.) 
Following Bcrghold et aZ4£ we define 

zi, n = (W n \e iGl - r \W n ), (A.3) 

and rewrite Eq. (|A.2|) as 

Nocc 6 
^ ' n=l 1=1 

where Gi = (1, 0, 0), G 2 = (0, 1, 0), G 3 = (0, 0, 1), G 4 = 
(1, 1, 0), G 5 = (1, 0, 1), and G 6 = (0, 1, 1), with respect 



The pseudo-wave-functions ty n are the variational quan- 
tities of the PAW method, and are expanded in recipro- 
cal space using plane waves. The index k is & shorthand 
for the atomic site R&, the angular momentum quantum 
numbers and m&, and the reference energy e^. The 
all-electron partial waves <pk are solutions of the radial 
Schrodinger equation for a non-spin-polarized reference 
atom at specific energies e& and for a specific angular 
momentum Ik- The pseudo-partial- waves cj>k are equiva- 
lent to the all-electron partial waves outside a core radius 
r c and match continuously onto <fik inside the core radius. 
The projector functions pk are dual to the pseudo-partial- 
waves (i.e., l$k\<t>l) = hi)- 

Following Ferretti et al.p^ we use Eqs. (|A.ip and (|A.5|) . 
to rewrite Eq. (|A.3|) as 



zi, n = {W n \e iG '- r \W n ) + J2(W n \pk){pi\W n ) 



x / e^' r fc (r)0,(r)-fo(r)<fc(r) dr 



(A.6) 



In our implementation, the integral over the partial waves 
in Eq. (|A.6[) is approximated by 



L 



&(r)e iG '- r dr, 



(A.7) 



where Qki( r ) is given by Eq. (27) of Ref. l5Cj, and Or,. 
is the PAW sphere at atomic site R&. [Note that the 
double sum in Eq. (|A.6|) is limited to be site-diagonal, 
i.e., Rfe = R;.] An exact expression for the integral in 
Eq. can be found in App. B of Ref. [zl 

Thus defined within the PAW formalism, the spread 
functional S[f7] is minimized with respect to U by means 
of the two-by-two orbital rotation method introduced by 
Edminston and Ruedenberg (see Sec. Ill A. of Ref. |46| 
and references therein). Once the maximally localized 
Wannier functions are constructed, we use 



^ 2 >» = 7il> - \zLn\ 2 ) + (v)l + 0(G% (A. 
^ ' 1=1 



and 



1 3 

(r)„ = (W n \v\W n ) = — J2 a i^ z i,n + 0(Gj), (A.9) 



i=i 
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where aj are the real space lattice vectors, to re- 
place the electronic part of the integrals J p m (r)rdr and 
/ p m (r)r 2 dr, in Eqs. Q and ©, by the sums J2i e ( r )i 
and J2i e ( r2 )i taken over the set of Wannier functions 
associated with molecular center m. 



Note that Eqs. (|A.8|) and (|A.9p are exact to order 
0(Gj), which for simple cubic supercells of lattice pa- 
rameter L, means the expressions for (r 2 ) n and (r)„ con- 
verge as L~ 2 in the limit of large L (see Fig. [3]). 
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